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This  paper  studies  the  application  of  discrete  multi-model  predictive  control  as  a  trajectory 
tracking  guidance  law  for  a  space  launcher.  Two  different  algorithms  are  developed,  each 
one  based  on  a  different  representation  of  launcher  translation  dynamics.  These  represen¬ 
tations  are  based  on  an  interpolation  of  the  linear  approximation  of  nonlinear  pseudo-five 
degrees  of  freedom  equations  of  translation  around  an  elliptical  Earth.  The  interpolation 
gives  a  linear-time-varying  representation  and  a  linear-fractional  representation.  They  are 
used  as  the  predictive  model  of  multi-model  predictive  controllers.  The  controlled  variables 
are  the  orbital  parameters,  and  constraints  on  a  terminal  region  for  the  minimal  accepted 
precision  are  also  included.  Use  of  orbital  parameters  as  the  controlled  variables  allows  for 
a  partial  definition  of  the  trajectory.  Constraints  can  also  be  included  in  multi-model 
predictive  control  to  reduce  the  number  of  unknowns  of  the  problem  by  defining  input 
shaping  constraints.  The  guidance  algorithms  are  tested  in  nominal  conditions  and  off- 
nominal  conditions  with  uncertainties  on  the  thrust.  The  results  are  compared  to  those  of  a 
similar  formulation  with  a  nonlinear  model  predictive  controller  and  to  a  guidance  method 
based  on  the  resolution  of  a  simplified  version  of  the  two-point  boundary  value  problem.  In 
nominal  conditions,  the  model  predictive  controllers  are  more  precise  and  produce  a  more 
optimal  trajectory  but  are  longer  to  compute  than  the  two-point  boundary  solution. 
Moreover,  in  presence  of  uncertainties,  developed  algorithms  exhibit  poor  robustness 
properties.  The  multi-model  predictive  control  algorithms  do  not  reach  the  desired  orbit 
while  the  nonlinear  model  predictive  control  algorithm  still  converges  but  produces  larger 
maneuvers  than  the  other  method. 

©  2013  IAA.  Published  by  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

It  is  well  known  that  there  are  two  approaches  to  guide 
a  vehicle  toward  its  final  destination  [1]:  predictor/correc¬ 
tor  methods  and  path  reference  methods.  Launch  ascent 
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guidance  of  a  space  launcher  is  no  different  than  any  other 
vehicle,  and  both  approaches  have  been  applied  to  it  in  the 
past  [2],  The  predictor/corrector  approaches  consist  in 
the  generation  of  a  new  trajectory  and  the  corresponding 
steering  commands  at  each  iteration  where  the  current  state 
is  the  initial  state  of  a  two-point  boundary  value  problem. 
Hence,  complex  optimization  algorithms  are  required  to 
solve  the  problem,  or  a  hypothesis  can  be  formulated  to 
simplify  the  problem.  For  the  launch  ascent  trajectory,  two 
hypotheses  on  the  Earth  gravity  approximation,  propor¬ 
tional  to  the  radius  [3]  or  uniform  [4],  give  an  analytical 
solution  of  the  costate  system  and  eliminate  the  need  for 
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complex  algorithms.  Both  of  these  hypotheses  are  well 
studied  and  can  be  considered  mature  techniques  suitable 
for  many  launch  applications  [4],  The  path  reference  meth¬ 
ods  compute  the  steering  commands  needed  to  follow  a 
predefined  trajectory.  As  presented  by  Shrivastava  et  al.  [2], 
many  linear  control  methods  with  time  scheduled  values 
have  been  studied  for  launch  ascent  guidance.  Among  them, 
Q-guidance  is  the  most  important  one  [5],  Recent  develop¬ 
ments  in  nonlinear  control  theory  boost  the  path  reference 
approaches  as  linearization  of  the  nonlinear  equations  of 
motion  is  not  required.  Nonlinear  model  predictive  control 
(NMPC)  [6],  neural  networks  [7]  and  sliding  mode  control 
[8]  are  examples  of  the  application  of  nonlinear  control 
methods  to  the  guidance  of  a  space  launcher. 

Even  if,  under  nominal  conditions,  the  trajectory  produced 
by  the  path  reference  methods  is  closer  to  the  optimal 
trajectory  than  the  trajectory  obtained  by  the  predictor/ 
corrector  methods  [2],  launch  ascent  guidance  is  traditionally 
achieved  by  predictor/corrector  methods  [9],  The  main  rea¬ 
son  for  this  choice  is  the  robustness  of  the  predictor/corrector 
methods  under  non-nominal  conditions;  however,  in  the 
exo-atmospheric  portion,  space  launchers  are  expected  to 
follow  their  nominal  trajectory  quite  closely  as  they  evolve  in 
an  environment  that  is  practically  perturbation-free.  There¬ 
fore,  path  reference  must  be  considered  for  this  portion  of  the 
launch,  mainly  for  vehicles  with  well-defined  missions  and 
reliable  components  [2], 

This  paper  focuses  on  MPC-based  path  reference 
approaches.  Through  integral  resolution,  continuous  time 
predictive  control  of  Lu  [6]  restricts  tracking  to  the 
in-plane  ascent  trajectory  around  a  spherical  Earth.  The 
proposed  guidance  laws  are  discrete  versions  of  predictive 
control.  In  opposition  to  continuous  predictive  control, 
discrete  formulation  has  a  finite  number  of  unknowns  and 
does  not  require  the  model  to  have  a  special  form  to  obtain 
the  solution.  Therefore,  the  in-plane  and  out-of-plane  ascent 
problem  of  a  motion  around  an  elliptical  Earth  can  be 
solved.  The  discrete  model  is  an  Euler  approximation  of 
the  pseudo-5  degrees  of  freedom  (pseudo-5DoF)  equations 
of  motion.  The  prediction  of  a  NMPC  is  made  directly  with 
these  equations  of  motion  while  the  varying  linear  repre¬ 
sentations  are  the  bases  of  multi-model  predictive  control 
(MMPC)  algorithms.  Vachon  et  al.  [10]  apply  discrete  NMPC 
to  space  launcher  exo-atmospheric  guidance.  The  varying 
linear  representations  are  time  interpolation  (linear-time- 
varying  representation  (ETVR)  and  linear-fractional  repre¬ 
sentation  (LFR))  of  a  set  of  linear  models  of  the  launcher 
translation  dynamics  [11], 

MMPC  algorithms  based  on  LTVR  and  LFR  are  well- 
known  approaches  used  to  handle  multiple  operation 
regimes  [12-15],  That  being  said,  they  are  mostly  imple¬ 
mented  as  linear  MPC  where  the  predictive  model  varies  at 
each  iteration  [12]  or  where  multiple  linear  MPC  algo¬ 
rithms  are  simultaneously  solved  and  a  second  algorithm 
weighs  the  results  based  on  the  current  operating  regime 
[14],  These  formulations  are  valid  for  systems  operating  in 
multiple  operating  regimes  but  where  the  variations  are 
slow  and  hence  the  prediction  around  a  single  operating 
regime  is  valid.  When  the  operating  regime  varies  inside 
the  prediction  horizon,  the  algorithms  are  based  on  min- 
max  optimization  to  obtain  the  varying  model  dynamically 


[13,15].  This  is  a  good  approach  when  the  sequence  of  the 
operating  regimes  is  not  known  a  priori.  In  a  launcher 
guidance  law,  the  prediction  must  be  made  over  multiple 
operating  regimes  that  are  known  a  priori.  Hence,  identify¬ 
ing  the  varying  linear  model  using  minmax  optimization  at 
each  time  step  is  not  necessary,  and  pre-identified  models 
can  be  used.  The  minmax  algorithm  is  therefore  converted 
into  a  single  minimization. 

Section  2  gives  an  overview  of  the  representations  of  the 
launcher  translation  dynamics  obtained  by  Vachon  et  al.  [11  ]. 
Section  3  introduces  the  NMPC  and  the  two  MMPC  formula¬ 
tions.  The  resulting  formulations  are  then  implemented  and 
tested  in  simulations.  Section  4  compares  the  results  with  a 
specific  predictor/corrector  method  [3]. 


2.  Equations  of  motion 

2.1.  Pseudo-five  degrees  of  freedom  equations 

As  stated  by  Zipfel  [16],  a  guidance  law  developed  for  a 
pseudo-5DoF  can  be  implemented  in  a  full  six-degree 
simulator  with  no  modifications.  Pseudo-5DoF  equations 
are  composed  of  the  three  degrees  of  freedom  equations 
of  translation  and  the  pseudo-two  degrees  of  freedom 
approximating  the  rotational  dynamics.  The  polar  coordi¬ 
nates  version  of  the  translation  equations  is  used  and  the 
controlled  rotational  dynamics  are  approximated  by  two 
first-order  transfer  functions  (one  for  the  in-plane  motion 
and  another  for  the  out-of-plane  motion).  A  time  constant 
of  1  s  for  both  transfer  functions  is  coherent  with  the 
requirements  of  an  exo-atmospheric  control  function  [17], 
Combining  these  two  sets  of  equations  gives  the  complete 
nonlinear  pseudo-5DoF  equations: 
m  =  Am  (la) 

f  =  v  sin (y)  (lb) 


= - m - &sin  7+S«  cos  x  cos  r 

-co2er  cos  <5(sin  5  cos  *  cos  y-  cos  5  sin  y) 
V  COS  y  cos  X 

v  cos  y  sin  / 
r  cos  <5 

T  cos  S  sin  q>  v  . 

= - -  +  -  sm  x  cos  y  tan  5 

mv  cos  y  r 

-  gs  sm  x  H - j — r  cos  s  sin  s  sin  x 

V  COS  y  V  COS  y 

+  -  cos  y+2coe  sin  x  cos  8 
coir  cos  5 

4-  — - (cos  /  cos  5 +  sin  S  cos  /  sm  y) 


(lc) 

(id) 

(le) 


(If) 


(lg) 
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S=Scom-S  (lh) 

<P  =  <Pcom-  <P  (1>) 

where  m  is  the  launcher  mass.  The  instantaneous  mass 
consumption  (Am)  and  the  thrust  magnitude  (T)  are  motor- 
related  constants.  The  r  and  v  variables  are  respectively  the 
magnitude  of  the  radius  (distance  between  the  Earth  center 
and  launcher  center  of  mass)  and  speed  vectors.  The  latitude 
(5),  the  longitude  (2),  the  flight-path  heading  (*•)  and  the 
flight-path  inclination  (/)  are  angles  describing  the  orienta¬ 
tion  of  the  radius  and  speed  vectors.  The  radius  vector  and  its 
corresponding  angles  are  presented  in  Fig.  1  while  the  speed 
vector  is  presented  in  Fig.  2.  Angles  8  and  q>  are  the  in-plane 
angle  and  out-of-plane  thrust  orientation  angles  respectively 
(Fig.  2).  This  representation  has  discontinuities  over  the  poles 
(<5=  +  90°)  but  they  are  not  harmful  to  the  sun-synchronous 
orbit  studied  later  in  this  work  (Section  4). 

In  these  equations,  the  variables  gr  and  gs,  representing 
the  gravitational  acceleration,  set  the  model  of  the  Earth. 
Most  of  space  launcher  guidance  functions  are  designed 
with  either  a  uniform  gravity  over  a  flat  Earth  or  a  gravity 
proportional  to  the  inverse  of  the  square  of  the  radius 
over  a  spherical  Earth  [2],  That  being  said,  launches  for 
high  inclination  orbits,  like  the  orbit  studied  as  part  of  this 
work,  are  expected  to  reach  high  latitudes  where  the 
Earth’s  ellipticity  has  a  significant  impact  on  the  gravity 
[18],  Its  inclusion  in  the  guidance  model  can  improve  the 
precision  of  the  resulting  trajectory  [19],  Hence,  the  second 


zonal  harmonic  is  considered  in  this  work.  The  corre- 


sponding  gravity  approximation  is 

gr=rl(l  +  1-5^(>)2(l-3sin25)) 

(2) 

sin  <5  cos  <5^ 

(3) 

Eq.  (1)  are  the  state  equations  of  the  nonlinear  model. 
A  second  set  of  equations,  the  output  equations,  is  needed 
to  complete  it.  Instead  of  tracking  the  speed  and  radius 
vectors  as  path  reference  methods  usually  do,  this  work 
suggests  tracking  the  orbital  parameters.  This  choice  is 
mainly  justified  by  the  possibility  of  partially  defining  the 
trajectory.  The  injection  anomaly,  namely  the  position  of 
the  launcher  on  its  orbit  during  the  injection,  is  irrelevant 
for  a  launcher.  In  this  work,  the  desired  orbit  is  circular  so 
the  argument  of  perigee  is  undefined  and  is  unnecessary 
to  the  trajectory  definition.  Moreover,  the  right  ascension 
of  the  ascending  node  is  also  omitted  but,  if  necessary, 
it  could  easily  be  added  to  the  output  vector.  Only  three 
parameters  of  the  six  related  to  the  components  of  the 
radius  and  speed  vectors  are  then  studied.  This  would 
not  be  possible  for  a  vector  tracking  algorithm  as  the 
six  components  are  needed  to  define  the  semi-major  axis, 
therefore,  fixing  the  five  other  orbital  parameters.  That 
being  said,  other  values  of  the  vectors  components  can 
give  the  desired  semi-major  axis  but  a  different  value  of  an 
unconsidered  orbital  parameter.  The  size  of  the  orbit  is 
normally  defined  by  geometric  constants:  semi-major  axis 
(a)  and  eccentricity  (e): 


6  V^r)  cos2  cos2  (5) 

The  subscript  /  refers  to  the  inertial  values  of  the  para¬ 
meters  based  on  the  inertial  speed  vector  (V/  =  v+me  x  r). 
Inertial  values  are  needed  because  the  orbital  parameters 
are  defined  in  an  inertial  frame  whereas  the  speed  vector 
of  Eq.  (1)  is  the  non-inertial  flight-path  speed.  That  being 
said,  numerical  problems  resulting  from  an  eccentricity 
close  to  zero,  arise  in  the  computation  of  the  linear  varying 
models  and  in  MPC  optimization  algorithms.  Therefore,  this 
work  suggests  defining  the  orbit  with  its  dynamic  constants: 
the  specific  angular  momentum  (h)  and  the  specific  mechan¬ 
ical  energy  (£): 

h  =  rvj  cos  r,  (6a) 

Z=2~y  (6b) 

The  derivatives,  inside  the  time-varying  representations, 
of  both  of  these  expressions  exist  and  give  finite  values. 
Furthermore,  both  values  are  far  from  being  null,  and  this 
will  help  in  the  MPC  algorithms. 

The  orbital  inclination  is  defined  by  the  angle  between 
the  orbital  plane  and  equatorial  plane: 
i  =  arccosf  sin  x\  cos  <5) 


(6c) 
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Eqs.  (1)  and  (6)  are  the  NMPC  predictive  model  and  the 
bases  of  the  linear  time  varying  representations.  These 
time-varying  representations  are  interpolations  of  a  family 
of  linear  models  obtained  from  a  referential  trajectory. 
Each  linear  model  is  an  evaluation  of  the  linear  tangent 
model.  These  representations,  LTVR  and  LFR,  are  presented 
in  the  next  two  sections. 

2.2.  Linear-time-varying  representation 

Time  is  the  scheduling  parameter  of  the  linear-time- 
varying  representation.  For  launcher  guidance,  using  time 
as  the  scheduling  parameter  is  not  problematic  as  the 
launcher  is  expected  to  follow  a  predefined  trajectory 
closely  where  time  is  a  meaningful  parameter.  Time  affects 
the  launch  trajectory  because  orbital  parameters  are 
linked  directly  to  it.  Furthermore,  for  the  application  in  a 
MMPC  algorithm  (Section  3),  time  must  be  measured 
or  estimated  as  it  defines  the  referential  trajectory.  Time 
evolution  is  also  known  a  priori,  which  simplifies  and 
accelerates  optimizations. 

LTVR  consists  in  representing  the  nonlinear  process 
(Eqs.  (1)  and  (6))  as  a  model  with  the  following  form: 
A„X=Ae(t)A0X+Be(t)A0ll  (7) 

A0y  =  Ce(t)A0x  (8) 

where  the  state  vector  (x)  is  [m  r  v  8  A  x  y  {>  q>]T,  the  input 
vector  (u)  is  the  commanded  thrust  orientation  angles 
([^com  <PcomiT )  and  the  output  vector  (y)  contains  the  studied 
orbital  parameters  ([h  £  i]T).  The  notation  A„  stands  for 
variations  around  the  operating  points  A0u  =  u-uop(t). 

The  matrices  Ae(t),  Be(t)  and  Q.(t)  are  the  interpolations 
of  a  family  of  linear  models.  They  consist  of  polynomial 
expressions  that  best  suit  the  family  of  linear  models. 
Given  the  large  differences  between  the  values  inside  the 
matrices,  polynomial  fitting  is  completed  separately  on 
each  element  of  the  matrices.  For  the  current  application, 
a  fifth  order  polynomial  is  suitable.  Hence,  for  example, 
the  polynomial  expression  of  the  element  located  in  the 
first  column  and  first  row  of  matrix  Ae  (Aell)  is  obtained 
by  solving  the  following  expression  of  a  least  squares 
problem: 


Aen  ' 

'  t<l)5  f(l)4  t(1)3  t®2  t(1)  1 

Aeu 

'A®' 

t(2)5  f<2)4  t(2)3  tpy2  t(2)  -] 

Aeh 

A® 

Aen 

jimf  t(m)4  t(m)3  t(m?  t(m)  } 

Ae|i 

A?f 

Ae°n 

(9) 

where  t(n)  is  the  value  of  t  at  the  nth  operating  point.  This 
equation  is  written  and  solved,  using  the  pseudo-inverse 
[20],  for  each  element  of  matrices  A  and  C.  Matrix  B  is 
constant. 

The  matrix  Ae(t)  in  the  Eq.  (7)  is  then  obtained  by 
evaluating  the  resulting  polynomial  expressions: 

Ae(t)=A°  +  A’t+-  +  A;;t5  (10) 


where  the  matrix  A°  is 


Ae°u  Aenu  -  Ae% 


Ae%x  Ae°g2  •••  Ae°g 


2.3.  Linear-fractional-representation 

Linear-fractional-representation  [21]  is  a  special  repre¬ 
sentation  of  a  linear-parameter-varying  (LPV)  system 
where  parameter  dependency  is  encapsulated  in  the  A 
block  of  the  linear-fractional  transformation  (LFT)  [22], 
Hence,  the  process  of  obtaining  a  LFR  is  similar  to  the 
previous  section.  Matrices  Ae(t),  Be(t)  and  C e(t)  are 
obtained  by  fitting  the  family  of  linear  models.  Then,  using 
the  generalized  Morton's  realization  technique  [23],  Eq.  (9) 
is  converted  in  the  M-A  form.  Polynomial  coefficients 
being  rational  functions  of  the  scheduling  parameters, 
the  generalized  Morton's  technique  is  needed.  Hence, 
having  the  minimal  order  LFR  is  no  longer  a  guarantee. 
The  matrices  of  Eqs.  (7)  and  (8)  are  obtained  by  computing 
the  resulting  upper-LFT  [22]: 

[  Ae(t)  Be(t)l 

L(t)  „  =M22  +  tM21I(I-tM11I)-,M21  (12) 

LFR  matrices  Mn,  M12,  M21  and  M22  result  from  a  least 
squares  problem.  Identity  matrix  (I)  has  the  same  dimen¬ 
sions  as  the  order  of  the  LFR.  Computing  the  LFR  requires 
inversing  a  matrix  in  which  dimension  is  the  same  as  the 
order  of  the  realization,  thus  the  importance  given  to  the 
order  of  the  representation.  The  two  representations  (LTVR 
and  LFR)  have  been  applied  and  validated  for  the  launch 
ascent  trajectory  [11], 

3.  Model  predictive  control 

Model  predictive  control  (MPC)  was  introduced  to  the 
control  theory  by  the  petro-chemical  industry  in  the 
late  1970s  [24],  It  is  a  mix  between  control  theory  and 
optimization  where  an  optimization  problem  is  defined 
based  on  the  control  requirements  to  be  met.  A  model  of 
the  system  calculates  the  predicted  future  outputs  (y)  in 
response  to  a  series  of  inputs  (u).  Based  on  the  minimiza¬ 
tion  of  a  cost  function  that  includes  the  difference  between 
the  desired  outputs  (d)  and  these  predicted  future  outputs, 
the  optimization  algorithm  finds  the  optimal  inputs  to 
be  applied  to  the  process.  Depending  on  the  predicting 
equations,  the  resulting  MPC  controller  can  be  a  linear 
MPC,  a  NMPC  or  MMPC.  The  optimal  values  are  then  sent 
to  the  system  and,  at  the  next  time  step,  the  optimization 
starts  over  using  the  new  measurements.  The  prediction 
horizon  then  keeps  being  shifted  forward.  This  is  the 
receding  horizon  principle.  Fig.  3  shows  an  example  of  a 
discrete  MPC  applied  on  a  continuous  linear  system. 

As  MPC  formulation  results  in  an  optimization  problem, 
constraints  can  be  included  in  the  definition  of  the  con¬ 
troller.  With  constraints,  physical  limitations  of  the  vehicle 
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Fig.  3.  Schematics  of  receding  horizon  principle. 


can  be  addressed  and  so  it  can  then  be  operated  closer  to 
its  optimal  performance  [25],  These  constraints  can  also 
define  a  terminal  region  ensuring  MPC  stability  [26]  or  be 
used  to  define  a  shape  for  the  input  vector  to  diminish 
computation  time  [27], 

That  being  said,  the  ability  to  handle  constraints  is  also 
a  drawback  of  MPC.  It  requires  the  resolution  of  a  con¬ 
strained  optimization  problem  where  convergence  inside 
limited  available  time  is  not  ensured,  leading  to  subopti- 
mal  or  diverging  solutions.  This  is  why  MPC  is  currently 
more  common  in  the  process  industry  where  settling 
times  are  longer  [25]  and  the  optimization  algorithm  has 
more  time  to  converge. 


3.1.  Formulation  for  launch  ascent  trajectory  tracking 

The  current  section  formulates  the  MPC  algorithm  for 
the  launch  ascent  guidance  problem.  The  MPC  criterion 
is  the  differences  between  the  desired  and  predicted 
outputs.  These  outputs,  the  orbital  parameters,  define  the 
trajectory,  and  the  MPC  becomes  a  path  reference  method 
where  the  launcher  tries  to  follow  a  predefined  trajectory. 
A  term  to  limit  drastic  changes  between  consecutive 
inputs  is  also  introduced  in  the  criterion: 

J=  i^dk+,-yk+i)T<l(dk+,-yk+l) 

+  ZAul-i  (13) 

where  the  predicted  model  outputs  (the  controlled  vari¬ 
ables)  are  the  orbital  parameters:  y  =  [h  £  i]1  and  the 
model  inputs  (the  manipulated  variables)  are  the  thrust 
orientation  angles  u=  [t9com  <pcom]r.  Vector  A  u  is  the  inputs 
increment  vector,  Auk+!  =  uk+ 1  —  uk+i  , .  MPC  produces  a 
combination  of  inputs  separated  in  time  that  results  in 
motions  that  are  instantaneously  unfeasible  and  makes  the 
control  of  underactuated  or  non-holonomic  systems  pos¬ 
sible  with  MPC  controllers  [28],  Prediction  and  control 
horizons,  respectively,  hp  and  hc,  and  the  weighting 


matrices  Q.  and  R  are  design  parameters  of  the  guidance 
law.  The  value  of  each  parameter  depends  on  the  studied 
trajectory  and  is  explained  later  in  Section  4. 

The  MPC  formulation  also  includes  constraints.  With  a 
prediction  horizon  long  enough  to  always  include  injection 
time,  terminal  orbital  parameters  are  constants  and 
can  define  constraints  ensuring  that  the  orbit  reached  is 
the  desired  one.  Furthermore,  using  inequality  constraints 
combined  with  a  precision  vector  ( ek )  allows  for  the 
minimal  accepted  precision  on  the  orbit  to  be  explicitly 
specified  in  the  guidance  law  initialization.  This  precision 
is  defined  in  terms  of  altitude  and  inclination.  Therefore, 
the  apsides  radius  (ra  and  rp)  are  more  relevant  than  the 
dynamics  constants  (h  and  £)  for  constraints  defining  orbit 


From  the  dynamics  constants,  apsides  radius 
as  follows: 


y/-  2c(l  —  h2) 


(14) 

1  calculated 


(15) 


-Me 

2e 


yj-2e(l-h2) 


(16) 


Input  shaping  constraints  are  also  implemented  by 
constraining  Auk+t  at  some  specific  index  l.  These  con¬ 
straints  give  a  predetermined  shape  to  the  input  vector 
resulting  in  a  reduction  of  the  number  of  unknowns  of  the 
optimization  algorithm  and  its  computational  time.  Again, 
these  are  design  parameters  and  their  values  will  be 
discussed  in  Section  4. 

Other  important  topics  to  study  are  the  availability  of 
the  state  vector  components  and  the  ability  of  the  gui¬ 
dance  law  to  eliminate  steady-state  errors.  In  a  complete 
guidance-navigation  and  control  scheme,  both  topics  are 
handled  by  the  navigation  function.  As  stated  by  Shi  et  al. 
[29],  the  navigation  function  acts  as  an  observer  that 
eliminates  bias  from  the  position  and  speed  vectors  and 
hence  eliminates  steady-state  errors.  Furthermore,  with 
this  observer,  all  components  become  estimated  and 
available  to  the  guidance  law. 


3.2.  Linear  multi-model  predictive  control 

The  last  element  defining  a  MPC  formulation  is  its 
predictive  model,  which  predicts  the  future  behaviour  of 
outputs.  This  model  includes  state  equations  and  output 
equations.  In  the  case  of  a  launch  ascent  trajectory  track¬ 
ing,  state  equations  are  Eq.  (1)  while  (2)  stands  for  the 
model  outputs.  As  those  equations  are  nonlinear,  the 
previous  formulation  is  a  NMPC.  This  NMPC  formulation 
was  successfully  applied  to  launch  ascent  trajectory  using 
these  nonlinear  equations  [10].  Flowever,  an  extension  to 
simplify  the  optimization  problem  and  accelerate  its 
resolution  may  be  an  interesting  alternative.  An  easy  way 
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to  do  it  is  to  use  a  linear  model  instead  of  a  nonlinear 
model.  That  being  said,  the  launcher's  operating  range 
makes  the  single  linear  model  unreliable  for  the  whole 
duration  of  the  launch,  but  the  variation  can  be  expressed 
as  a  LTVR  [11],  The  idea  of  the  proposed  formulation  is  to 
modify  the  NMPC  by  converting  the  nonlinear  equations  to 
LTVR  (Section  2.2)  or  LFR  (Section  2.3)  in  order  to  obtain  a 
MMPC  with  pre-identified  models.  Furthermore,  as  the 
models  are  pre-identified  and  the  scheduling  parameters 
are  known  a  priori  the  minmax  optimization  of  a  MMPC 
predicting  around  multiple  operating  regimes  like  Lu  and 
Arkun  [13]  is  converted  into  a  single  minimization.  There¬ 
fore,  nonlinear  model  predictive  control  is  a  special  case  of 
MPC  in  which  the  predicting  equations  are  nonlinear  and 
MMPC  is  a  also  a  special  case  of  MPC  in  which  the  equations 
are  either  the  LTVR  or  the  LFR. 

3.3.  Stability  issues  of  MPC 

A  closed-loop  system  is  stable  at  a  set  point  if,  when 
starting  out  near  this  point,  the  system  stays  close  to  this 
set  point  forever  [30],  In  the  case  of  launcher  guidance,  the 
set  point  corresponds  to  the  desired  orbit,  and  the  stability 
of  the  system  at  this  orbit  ensures  its  sustainability.  MPC 
stability  resides  in  problem  feasibility  and  cost  function 
monotonicity  [30],  The  previous  formulation  of  the  pro¬ 
blem  is  somewhat  similar  to  the  MPC  dual  mode  control. 
This  formulation  includes  a  constant  ellipsoid  around  the 
set  point,  defined  by  constraints  on  yk+h  ,  to  define  two 
operation  modes.  Outside  this  ellipsoid,  an  MPC  controller 
is  employed  and,  once  inside  the  ellipsoid,  the  control  is 
switched  to  a  static  linear  state  feedback.  As  demonstrated 
by  Kwon  and  Han  [30]  for  linear  systems  and  by  Grime 
and  Pannek  [26]  for  nonlinear  cases,  this  dual-mode 
formulation  ensures  the  cost  monotonicity  condition,  de 
Nicolao  et  al.  [31  ]  proved  the  MMPC  stability  with  constant 
terminal  region.  In  the  current  formulation,  with  tk  in  the 
constraints,  the  terminal  region  varies  and  the  dynamics 
of  the  system  inside  also  vary.  That  being  said,  if  stability 
using  static  linear  state  feedback  is  possible  inside  the 
largest  region,  MMPC  stability  is  also  possible.  Dynamic 
variations  are  known  and  bounded;  it  is  then  possible 
to  verify  that  a  static  linear  feedback  is  robust  to  these 
variations  ensuring  worst-case  stability.  However,  these 
stability  proofs  are  based  on  the  hypothesis  that  the 
predictive  model  is  a  perfect  representation  of  the  system 
dynamics,  which  is  not  the  case  for  the  current  application. 

Furthermore,  even  if  cost  monotonicity  can  be  proven, 
the  finite  prediction,  control  horizons  and  additional  con¬ 
straints  makes  satisfaction  of  the  terminal  constraints 
more  difficult  to  guarantee  [28],  The  problem  feasibility 
condition,  the  second  requirement  to  guarantee  stability, 
cannot  therefore  be  guaranteed;  hence,  neither  can  the 
MPC  stability. 


4.  Simulations  and  results  analysis 

NMPC  and  MMPC  trajectory  tracking  laws  defined  in 
the  previous  section  are  applied  to  the  exo-atmospheric 
trajectory  of  a  three-stages-to-orbit  launcher  to  reach  a 


circular  sun-synchronous  orbit  with  a  semi-major  axis 
of  6871  km  and  an  orbital  inclination  of  97.375°.  Before 
applying  both  MPC  formulations,  the  referential  trajectory 
is  obtained  by  solving  an  off-line  trajectory  generation 
algorithm  [19],  The  trajectory  obtained  is  a  burn-coast- 
burn  one  where  the  coast  arc  occurs  immediately  after  the 
second-stage  jettison. 

With  the  separation  of  the  trajectory  by  the  coast 
phase,  the  size  of  the  NMPC  and  MMPC  formulations  can 
be  diminished.  During  the  coast  phase,  the  launcher  is 
unpropelled  and  the  orbit  cannot  be  modified.  Therefore, 
the  complete  trajectory  of  the  launcher  is  divided  into 
two  distinct  guided  phases  and  MPC  problems.  During  the 
second-stage  burn,  the  guidance  law  tries  to  inject  the 
third  stage  and  the  attached  payload  on  the  coast  orbit, 
calculated  by  the  off-line  algorithm.  When  the  third  stage 
is  active,  the  objective  of  the  guidance  law  is  to  put  the 
payload  on  the  desired  orbit  for  its  mission  (Fig.  4). 

With  this  division  into  two  consecutive  phases,  some  of 
the  parameters  differ  whether  the  second  or  third  stage  is 
active.  These  parameters  are  the  prediction  and  control 
horizons  and  the  terminal  region.  During  the  second-stage 
burn,  the  terminal  region  is  defined  around  the  coast  orbit 
and  around  the  desired  orbit  for  the  third-stage  burn. 
As  the  constraints  (Eq.  (14))  are  defined  at  the  end  of  the 
prediction  horizon,  the  latter  must  be  long  enough  to 
always  include  the  injection  time.  For  an  undivided  tra¬ 
jectory,  the  minimal  duration  of  the  prediction  would  be 
508.2  s,  while  for  the  divided  trajectory,  the  minimal 
durations  are  respectively  48.2  and  104.4  s  for  the  pre¬ 
coast  and  the  post-coast  problems.  In  MPC,  when  the 
control  horizon  is  shorter  than  the  prediction  horizon,  all 
input  increments  after  the  control  horizon  are  set  to  zero. 
That  being  said,  as  demonstrated  by  the  off-line  algorithm 
[19],  the  optimal  thrust  orientation  needed  to  reach  the 
orbit  changes  all  the  way  up  to  the  injection.  To  reproduce 
this  optimal  shape,  the  algorithm  must  therefore  be  able  to 
modify  the  thrust  orientation  for  the  whole  prediction 
horizon.  Hence,  the  control  horizon  must  be  as  long  as  the 
prediction  horizon.  A  shorter  control  horizon  means  fewer 
unknowns  for  the  optimization  problem  and  should  trans¬ 
late  into  a  shorter  computational  time.  For  the  current 
application,  with  a  1  s  guidance  period,  both  the  prediction 
and  control  horizons  are  49  for  the  pre-coast  burn  and  105 
steps  for  the  post-coast  phase  as  opposed  to  the  509  for 
the  undivided  trajectory.  The  pre-coast  burn  phase  starts 
with  the  fairing  jettison,  170.8  s  after  the  launch,  and  ends 
at  the  beginning  of  the  coast  phase  (219  s).  As  for  the  post¬ 
coast  phase,  it  starts  at  the  end  of  the  coast  phase  (569.4  s) 


Fig.  4.  Launch  trajectory  separation  in  two  consecutive  MPC  problems. 
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and  ends  with  the  injection  of  the  payload  on  its  final  orbit 
(673.9  s). 

As  the  order  of  magnitude  of  the  output  vector  is 
similar  in  both  phases,  weighing  matrices  can  be  the  same 
for  the  two  phases.  The  values  of  these  diagonal  matrices 
(Q.  and  R)  eliminate  problems  resulting  from  different 
orders  of  magnitude  inside  the  criterion.  The  implicit 
unit  of  each  element  gives  a  good  idea  of  the  order  of 
magnitude  of  the  matrix  components.  Hence,  the  element 
of  Q.  related  to  the  specific  angular  momentum  should  be 
low  to  bring  the  three  orbital  parameters  to  the  same 
order  of  magnitude.  These  matrices  set  the  relative  impor¬ 
tance  of  each  parameter.  Therefore,  as  they  are  more 
fuel  consuming,  out-of-plane  maneuvers  should  be  hardly 
penalized.  This  is  done  by  setting  a  bigger  weight  to  the 
out-of-plane  angles  and  the  orbital  inclination.  The  ele¬ 
ment  of  R  related  to  the  out-of-plane  angles  (R22)  and  the 
parameter  of  Q.  related  to  the  orbit  inclination  (Q33)  should 
be  higher  than  when  the  matrices  serve  only  for  normal¬ 
ization.  Preliminary  tests  showed  that  the  following  matrices 
Q  and  R  yield  good  results: 

["  1  x  10“6  0  0  1 

Q.Jj  0  1  X  102  0  (17) 

[  0  0  1  x  103J 


"fV05 

Concerning  the  constraints,  the  desired  orbital  para¬ 
meters  are  set  to  either  the  coast  orbital  parameters  or 
the  injection  orbital  parameters  depending  on  the  phase. 
Given  that  they  try  to  counteract  the  disturbances  in  a 
short  amount  of  time,  guidance  algorithms  which  do  not 
handle  thrust  cut-off  and  request  near-perfect  injection 
result  in  large  maneuvers  close  to  injection  [32],  For  the 
current  formulations,  this  problematic  behavior  can  be 
avoided  by  implementing  a  time-varying  precision  vector 
(e/t).  Earlier  in  the  launch,  it  is  set  at 

c0  =  [lxlO"'  1x10"’  lxl0“3]r  (19) 

to  request  a  near-perfect  launch  and,  as  the  remaining 
time  before  injection  decreases,  it  constantly  increases  to 
reach  the  desired  precision  for  each  orbital  parameter  at 
injection  time,  which,  for  the  current  application,  is  set  at 

ehp=[4  4  1  xlCrY  (20) 

The  implemented  time-varying  scheme  for  the  precision  is 
then 


ie0  if  k  <  ftp  -  25 

evor  if  hp-25  <k<hp-5  (21) 

ehp  if  fc  >  Zip  -  5 


of  the  horizon  followed  by  sequences  of  one  free  input 
increment  and  three  input  increments  constrained  to  be  null. 
In  the  pre-coast  phase,  the  tested  shape  is  composed  of  9 
free  increments  followed  by  10  sequences  (Fig.  5).  As  for  the 
shape  of  the  post-coast  phase,  it  consists  of  41  free  inputs 
and  16  sequences  (Fig.  6). 

Guidance  laws  are  applied  to  a  launch  simulator.  In  the 
first  test  (Section  4.1 ),  only  random  noises  corresponding 
to  the  output  precision  of  an  integrated  GPS/INS  naviga¬ 
tion  function  are  present.  Based  on  the  results  of  Maki 
[33],  random  white  noises  with  a  covariance  of  35  m  on 
the  position  and  0.1  m/s  on  the  velocity  are  reasonable 
estimates  of  an  achievable  precision  for  a  space  launcher 
navigation  system.  This  test  is  conducted  to  investigate 
guidance  law  performances  in  nominal  operating  condi¬ 
tions.  As  the  proposed  guidance  laws  are  path  reference 
methods,  they  should  perform  well  in  this  case  [2],  The 
second  series  of  tests  (Section  4.2)  add  motor  uncertainties 
to  the  random  noises  of  the  first  case.  The  uncertainties  are 
modeled  by  a  difference  between  the  thrust  magnitude  of 
the  simulator  and  the  one  of  the  predictive  model  of  the 
guidance  laws.  For  the  two  series  of  tests,  the  proposed 
MMPCs  using  the  LTVR  (LTVR-MPC)  and  the  LFR  (LFR- 
MPC),  are  compared  to  the  NMPC  algorithm  [10]  modified 
to  handle  the  time-varying  precision  and  to  a  specific 
predictor/corrector  method  [3], 


where  evar  is  the  equation  of  the  increasing  slope 


4.1.  Simulation  with  random  noises  only 


evar  =  (fc-  20 hp)  +  25ehp  -  5e0  (22) 

Input  shaping  constraints  also  give  a  predetermined 
shape  to  the  input  increments,  which  translates  into  a 
reduction  in  the  number  of  unknowns  in  the  optimization. 
This  shape  consists  of  free  input  increments  at  the  beginning 


Fig.  7  presents  the  applied  commanded  thrust  orienta¬ 
tions  obtained  with  the  four  algorithms  applied  to  a  simula¬ 
tion  with  only  random  noises  as  perturbations.  Tables  1 
and  2  present  the  orbits  reached  at  the  end  of  the  two  phases. 

As  expected,  in  nominal  operating  conditions,  path  refer¬ 
ence  algorithms  based  perform  better  than  the  comparison 
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Table  1 

Reached  coast  orbits  for  random  noises  only. 


Guidance  method  Apogee  Perigee  Orbital 

radius  (km)  radius  (km)  inclination  (deg) 


Desired  6862.5 

NMPC  6862.6 

LTVR-MPC  6862.5 

LFR-MPC  6862.5 

Predictor/corrector  6861.2 


1294.0  97.232 

1294.9  97.233 

1295.0  97.233 

1294.8  97.230 

1304.2  97.308 


Table  2 

Reached  final  orbits  for  random  noises  only. 


Guidance  method  Apogee  Perigee  Orbital 

radius  (km)  radius  (km)  inclination  (deg) 


Desired  6871.0  6871.0 

NMPC  6872.9  6869.9 

LTVR-MPC  6875.0  6867.7 

LFR-MPC  6872.4  6870.3 

Predictor/corrector  6871.3  6867.3 


97.375 

97.376 

97.375 

97.376 
97.375 


predictor/corrector  method.  Commanded  thrust  angles  are 
smaller  and  the  third  stage  duration  is  shorter.  This  translates 
into  a  maximal  payload  mass  of  187.5  kg  for  the  MPC 
algorithms  and  186.2  kg  for  the  comparison  predictor/correc¬ 
tor  method.  The  orbit  reached  by  the  three  MPC  algorithms 
are,  at  worst,  as  precise  as  the  referential  method.  There  are 
no  significant  differences  between  the  LFR,  LTVR  and  non¬ 
linear  model  [11];  this  is  also  visible  in  Fig.  7.  The  angles 
obtained  are  nearly  identical  for  the  three  algorithms.  That 
being  said,  the  investigation  of  computation  time  (Fig.  8) 
draws  a  completely  different  conclusion. 

Even  if  it  would  be  logical  to  think  that  path  reference 
algorithms  should  have  a  better  computation  time,  MPC 
formulations  are  worse  than  the  comparison  predictor/ 
corrector  method  by  far.  This  can  be  explained  by  the  MPC 
formulation.  It  results  in  a  constrained  optimization  that  is 
not  easier  to  solve  than  the  predictor/corrector  method. 
The  complexity  of  MPC  formulations  explains  why  the 


predictor/corrector  method,  which  uses  hypothesis  to 
simplify  the  optimization  problem,  is  quicker.  The  differ¬ 
ence  between  the  three  MPCs  is  explained  by  the  complex¬ 
ity  of  predictive  models.  The  LTVR  is  the  easiest  to 
compute,  the  LFR,  by  the  matrix  inversion,  is  longer  and 
the  NMPC  is  the  longest. 


4.2.  Simulation  with  random  noises  and  motors  misfunction 

In  this  section,  multiple  tests  were  conducted  to  verify 
the  robustness  of  the  developed  guidance  laws  to  a  model 
mismatch.  To  do  so,  the  thrust  of  the  real  launcher  was 
modified.  The  main  results  of  these  tests  are  the  ineffi¬ 
ciency  of  the  NMPC  and,  particularly  the  MMPCs.  When 
the  malfunction  is  relatively  small,  roughly  below  3%  for 
less  than  40  s,  the  four  algorithms  are  almost  unaffected  by 
the  model  error;  however,  when  the  difference  is  larger, 
the  MMPCs  have  trouble  providing  the  desired  orbit. 
Table  3  presents  the  final  orbit  achieved  when  the  third 
stage  thrust  is  3%  lower  than  its  nominal  value  for  its 
whole  duration  and  Table  4  shows  the  results  when  a  3% 
malfunction  occurs  640  s  after  the  launch  and  lasts  until 
the  injection. 

These  results  can  be  explained  by  the  fact  that  the 
linear  MMPC  approximations  are  valid  only  around 
the  referential  trajectory  for  which  they  have  been 
defined.  That  being  said,  a  malfunction  steers  them  away 
from  this  trajectory.  The  validity  of  the  approximation  is 
gradually  lost,  decreasing  performances.  As  for  the  NMPC, 
it  is  not  defined  around  a  single  trajectory  and  is  valid  for  a 
wider  range.  So  it  does  converge  to  the  desired  orbit; 
however,  it  requires  large  maneuvers  to  do  so  because  it 
tries  to  follow  a  trajectory  that  cannot  be  attained  with  the 
launcher's  new  configuration.  The  comparison  predictor/ 
corrector  method  seems  to  be  less  sensitive  to  motor 
malfunctions;  its  result  is  relatively  unaffected  by  a  change 
in  the  thrust  magnitude.  This  is  explained  by  the  fact  that 
the  design  of  this  method  is  based  on  a  highly  simplified 
model.  Hence,  the  algorithm  must  be  robust  to  model 
uncertainties  even  in  nominal  conditions.  These  results  are 
coherent  with  the  preference  for  the  predictor/corrector 
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Fig.  8.  Relative  computation  time  needed  for  random  noises  only. 


Table  3 

Reached  final  orbits  for  random  noises  and  a  3%  thrust  difference  for 
105  s. 


Guidance  method  Apogee  Perigee  Orbital 

radius  (km)  radius  (km)  inclination  (deg) 


Desired  6871.0 

NMPC  6868.0 

LTVR-MPC  6871.9 

LFR-MPC  6870.1 

Predictor/corrector  6871.1 


97.375 

97.369 

97.363 

97.358 

97.374 


representation,  which  is  significantly  quicker.  That  being  said, 
its  computation  time  still  underperforms  a  predictor/correc¬ 
tor  method  used  for  comparison.  Furthermore,  from  the 
simulation  results,  it  seems  that  to  become  a  truly  viable 
solution  for  a  space  launcher  exo-atmospheric  guidance 
law,  the  multi-model  predictive  control  algorithms  probably 
still  need  to  have  their  robustness  to  input  uncertainties 
improved. 


Acknowledgments 


Table  4 

Reached  final  orbits  for  random  noises  and  a  3%  thrust  difference  for 
34.9  s. 


Guidance  method  Apogee  Perigee  Orbital 

radius  (km)  radius  (km)  inclination  (deg) 


Desired  6871.0  6871.0 

NMPC  6872.2  6870.8 

LTVR-MPC  6874.4  6868.7 

LFR-MPC  6873.8  6868.5 

Predictor/corrector  6871.8  6869.8 


97.375 

97.376 
97.375 
97.375 
97.375 


methods  in  space  launcher  guidance  because  they  perform 
better  under  off-nominal  conditions  [9]. 

5.  Conclusion 

This  work  has  developed  trajectory  tracking  guidance 
laws  based  on  discrete  model  predictive  control.  A  first 
algorithm  is  based  on  a  nonlinear  model  of  the  translation 
dynamics,  and  the  others  use  a  time-varying  approximation 
of  the  nonlinear  translation  dynamics  as  their  predictive 
model.  Tracking  is  done  using  the  orbital  parameters  as 
outputs;  this  enables  a  partial  definition  of  the  orbit.  The 
orbital  parameters  define  an  invariant  terminal  zone  using 
constraints  on  the  last  prediction.  This  is  similar  to  the 
algorithm  in  the  dual  mode  MPC,  which  ensures  stability 
and  makes  it  possible  to  include  the  desired  precision  directly 
in  the  definition  of  the  guidance  law.  Moreover,  to  reduce 
the  number  of  unknowns  in  the  optimization  algorithm, 
the  bum-coast-burn  trajectory  is  divided  into  two  separate 
guidance  problems,  and  input  shaping  constraints  are  imple¬ 
mented.  Multi-model  predictive  control  algorithms  prove  to 
be  a  viable  alternative  to  nonlinear  model  predictive  control 
for  a  nominal  launch,  particularly  the  linear-time-varying 
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